!\Description:
!  Reverse communication interface for applying NP additional steps to
!  a K step nonsymmetric Arnoldi factorization.
!
!  Input:  OP*V_{k}  -  V_{k}*H = r_{k}*e_{k}^T
!
!          with (V_{k}^T)*B*V_{k} = I, (V_{k}^T)*B*r_{k} = 0.
!
!  Output: OP*V_{k+p}  -  V_{k+p}*H = r_{k+p}*e_{k+p}^T
!
!          with (V_{k+p}^T)*B*V_{k+p} = I, (V_{k+p}^T)*B*r_{k+p} = 0.
!
!  where OP and B are as in znaupd.  The B-norm of r_{k+p} is also
!  computed and returned.
!\Arguments
!  IDO     Integer.  (INPUT/OUTPUT)
!          Reverse communication flag.
!          -------------------------------------------------------------
!          IDO =  0: first call to the reverse communication interface
!          IDO = -1: compute  Y = OP * X  where
!                    IPNTR(1) is the pointer into WORK for X,
!                    IPNTR(2) is the pointer into WORK for Y.
!                    This is for the restart phase to force the new
!                    starting vector into the range of OP.
!          IDO =  1: compute  Y = OP * X  where
!                    IPNTR(1) is the pointer into WORK for X,
!                    IPNTR(2) is the pointer into WORK for Y,
!                    IPNTR(3) is the pointer into WORK for B * X.
!          IDO =  2: compute  Y = B * X  where
!                    IPNTR(1) is the pointer into WORK for X,
!                    IPNTR(2) is the pointer into WORK for Y.
!          IDO = 99: done
!          -------------------------------------------------------------
!          When the routine is used in the "shift-and-invert" mode, the
!          vector B * Q is already available and do not need to be
!          recomputed in forming OP * Q.
!
!  BMAT    Character*1.  (INPUT)
!          BMAT specifies the type of the matrix B that defines the
!          semi-inner product for the operator OP.  See znaupd.
!          B = 'I' -> standard eigenvalue problem A*x = lambda*x
!          B = 'G' -> generalized eigenvalue problem A*x = lambda*M**x
!
!  N       Integer.  (INPUT)
!          Dimension of the eigenproblem.
!
!  K       Integer.  (INPUT)
!          Current size of V and H.
!
!  NP      Integer.  (INPUT)
!          Number of additional Arnoldi steps to take.
!
!  NB      Integer.  (INPUT)
!          Blocksize to be used in the recurrence.
!          Only work for NB = 1 right now.  The goal is to have a
!          program that implement both the block and non-block method.
!
!  RESID   complex(8) array of length N.  (INPUT/OUTPUT)
!          On INPUT:  RESID contains the residual vector r_{k}.
!          On OUTPUT: RESID contains the residual vector r_{k+p}.
!
!  RNORM   real(8) scalar.  (INPUT/OUTPUT)
!          B-norm of the starting residual on input.
!          B-norm of the updated residual r_{k+p} on output.
!
!  V       complex(8) N by K+NP array.  (INPUT/OUTPUT)
!          On INPUT:  V contains the Arnoldi vectors in the first K
!          columns.
!          On OUTPUT: V contains the new NP Arnoldi vectors in the next
!          NP columns.  The first K columns are unchanged.
!
!  LDV     Integer.  (INPUT)
!          Leading dimension of V exactly as declared in the calling
!          program.
!
!  H       complex(8) (K+NP) by (K+NP) array.  (INPUT/OUTPUT)
!          H is used to store the generated upper Hessenberg matrix.
!
!  LDH     Integer.  (INPUT)
!          Leading dimension of H exactly as declared in the calling
!          program.
!
!  IPNTR   Integer array of length 3.  (OUTPUT)
!          Pointer to mark the starting locations in the WORK for
!          vectors used by the Arnoldi iteration.
!          -------------------------------------------------------------
!          IPNTR(1): pointer to the current operand vector X.
!          IPNTR(2): pointer to the current result vector Y.
!          IPNTR(3): pointer to the vector B * X when used in the
!                    shift-and-invert mode.  X is the current operand.
!          -------------------------------------------------------------
!
!  WORKD   complex(8) work array of length 3*N.  (REVERSE COMMUNICATION)
!          Distributed array to be used in the basic Arnoldi iteration
!          for reverse communication.  The calling program should not
!          use WORKD as temporary workspace during the iteration !!!!!!
!          On input, WORKD(1:N) = B*RESID and is used to save some
!          computation at the first step.
!
!  INFO    Integer.  (OUTPUT)
!          = 0: Normal exit.
!          > 0: Size of the spanning invariant subspace of OP found.
!\Remarks
!  The algorithm implemented is:
!
!  restart = .false.
!  Given V_{k} = [v_{1}, ..., v_{k}], r_{k};
!  r_{k} contains the initial residual vector even for k = 0;
!  Also assume that rnorm = || B*r_{k} || and B*r_{k} are already
!  computed by the calling program.
!
!  betaj = rnorm ; p_{k+1} = B*r_{k} ;
!  For  j = k+1, ..., k+np  Do
!     1) if ( betaj < tol ) stop or restart depending on j.
!        ( At present tol is zero )
!        if ( restart ) generate a new starting vector.
!     2) v_{j} = r(j-1)/betaj;  V_{j} = [V_{j-1}, v_{j}];
!        p_{j} = p_{j}/betaj
!     3) r_{j} = OP*v_{j} where OP is defined as in znaupd
!        For shift-invert mode p_{j} = B*v_{j} is already available.
!        wnorm = || OP*v_{j} ||
!     4) Compute the j-th step residual vector.
!        w_{j} =  V_{j}^T * B * OP * v_{j}
!        r_{j} =  OP*v_{j} - V_{j} * w_{j}
!        H(:,j) = w_{j};
!        H(j,j-1) = rnorm
!        rnorm = || r_(j) ||
!        If (rnorm > 0.717*wnorm) accept step and go back to 1)
!     5) Re-orthogonalization step:
!        s = V_{j}'*B*r_{j}
!        r_{j} = r_{j} - V_{j}*s;  rnorm1 = || r_{j} ||
!        alphaj = alphaj + s_{j};
!     6) Iterative refinement step:
!        If (rnorm1 > 0.717*rnorm) then
!           rnorm = rnorm1
!           accept step and go back to 1)
!        Else
!           rnorm = rnorm1
!           If this is the first time in step 6), go to 5)
!           Else r_{j} lies in the span of V_{j} numerically.
!              Set r_{j} = 0 and rnorm = 0; go to 1)
!        EndIf
!  End Do
!-----------------------------------------------------------------------
subroutine cnaitr(ido, n, k, np, resid, rnorm, v, h, ldh, ipntr, workd, info)
    implicit none
    integer ido, info, k, ldh, n, np, ipntr(3)
    complex h(ldh, np), resid(n), v(n, np), workd(3 * n)
    real rnorm
    ! Parameters
    complex, parameter :: one = (1.0, 0.0), zero = (0.0, 0.0)
    real, parameter :: rone = 1.0
    logical first, orth1, orth2, rstart, step3, step4
    integer ierr, i, infol, ipj, irj, ivj, iter, itry, j, jj
    real    ovfl, smlnum, tst1, ulp, unfl, betaj, temp1, rnorm1, wnorm
    save    first, orth1, orth2, rstart, step3, step4, ierr, ipj, irj, ivj, iter, itry, j, ovfl, &
            betaj, rnorm1, smlnum, ulp, unfl, wnorm
    ! External Functions
    complex, external :: cdotc
    real, external :: slamch, scnrm2, clanhs, slapy2
    ! Data statements
    data       first / .true. /

    if (first) then
        ! Set machine-dependent constants for the the splitting and deflation criterion.
        ! If norm(H) <= sqrt(OVFL), overflow should not occur.
        ! REFERENCE: LAPACK subroutine zlahqr
        unfl = slamch('safe minimum')
        ovfl = real(1.0 / unfl)
        call slabad(unfl, ovfl)
        ulp = slamch('precision')
        smlnum = unfl * (n / ulp)
        first = .false.
    end if

    if (ido == 0) then
        ! Initial call to this routine
        info = 0
        step3 = .false.
        step4 = .false.
        rstart = .false.
        orth1 = .false.
        orth2 = .false.
        j = k + 1
        ipj = 1
        irj = ipj + n
        ivj = irj + n
    end if
    ! When in reverse communication mode one of:
    ! STEP3, STEP4, ORTH1, ORTH2, RSTART will be .true. when ....
    ! STEP3: return from computing OP*v_{j}.
    ! STEP4: return from computing B-norm of OP*v_{j}
    ! ORTH1: return from computing B-norm of r_{j+1}
    ! ORTH2: return from computing B-norm of correction to the residual vector.
    ! RSTART: return from OP computations needed by zgetv0.
    if (step3)  go to 50
    if (step4)  go to 60
    if (orth1)  go to 70
    if (orth2)  go to 90
    if (rstart) go to 30
    ! Else this is the first step
    !--------------------------------------------------------------%
    !        A R N O L D I     I T E R A T I O N     L O O P       |
    ! Note:  B*r_{j-1} is already in WORKD(1:N)=WORKD(IPJ:IPJ+N-1) |
    !--------------------------------------------------------------%
1000 continue
    !---------------------------------------------------%
    ! STEP 1: Check if the B norm of j-th residual      |
    ! vector is zero. Equivalent to determine whether   |
    ! an exact j-step Arnoldi factorization is present. |
    !---------------------------------------------------%
    betaj = rnorm
    if (rnorm > 0.0) go to 40
    ! ITRY is the loop variable that controls the maximum amount of times that a restart is attempted.
    betaj = 0.0
    itry = 1
20  continue
    rstart = .true.
    ido = 0
30  continue
    ! If in reverse communication mode and RSTART = .true. flow returns here.
    call cgetv0(ido, .false., n, j, v, resid, rnorm, workd, ierr)
    if (ido /= 99) return
    if (ierr < 0) then
        itry = itry + 1
        if (itry <= 3) go to 20
        ! Give up after several restart attempts.
        ! Set INFO to the size of the invariant subspace  which spans OP and exit.
        info = j - 1
        ido = 99
        return
    end if
40  continue
    !---------------------------------------------------------%
    ! STEP 2:  v_{j} = r_{j-1}/rnorm and p_{j} = p_{j}/rnorm  |
    ! Note that p_{j} = B*r_{j-1}. In order to avoid overflow |
    ! when reciprocating a small RNORM, test against lower    |
    ! machine bound.                                          |
    !---------------------------------------------------------%
    call ccopy(n, resid, 1, v(1, j), 1)
    if (rnorm >= unfl) then
        temp1 = 1.0 / rnorm
        call csscal(n, temp1, v(1, j), 1)
        call csscal(n, temp1, workd(ipj), 1)
    else
        ! To scale both v_{j} and p_{j} carefully use LAPACK routine zlascl
        call clascl('General', i, i, rnorm, rone, n, 1, v(1, j), n, infol)
        call clascl('General', i, i, rnorm, rone, n, 1, workd(ipj), n, infol)
    end if
    ! STEP 3:  r_{j} = OP*v_{j}; Note that p_{j} = B*v_{j}
    ! Note that this is not quite yet r_{j}. See STEP 4
    step3 = .true.
    call ccopy (n, v(1, j), 1, workd(ivj), 1)
    ipntr(1) = ivj
    ipntr(2) = irj
    ipntr(3) = ipj
    ido = 1
    ! Exit in order to compute OP*v_{j}
    return
50  continue
    ! Back from reverse communication; WORKD(IRJ:IRJ+N-1) := OP*v_{j} if step3 = .true.
    step3 = .false.
    ! Put another copy of OP*v_{j} into RESID.
    call ccopy(n, workd(irj), 1, resid, 1)
    ! STEP 4:  Finish extending the Arnoldi factorization to length j.
    call ccopy(n, resid, 1, workd(ipj), 1)
60  continue
    ! Back from reverse communication;
    ! WORKD(IPJ:IPJ+N-1) := B*OP*v_{j}
    ! if step4 = .true.
    step4 = .false.
    ! The following is needed for STEP 5. Compute the B-norm of OP*v_{j}.
    wnorm = scnrm2(n, resid, 1)
    ! Compute the j-th residual corresponding
    ! to the j step factorization.
    ! Use Classical Gram Schmidt and compute:
    ! w_{j} <-  V_{j}^T * B * OP * v_{j}
    ! r_{j} <-  OP*v_{j} - V_{j} * w_{j}

    ! Compute the j Fourier coefficients w_{j} WORKD(IPJ:IPJ+N-1) contains B*OP*v_{j}.
    call cgemv('C', n, j,  one, v, n, workd(ipj), 1, zero, h(1, j), 1)
    !Orthogonalize r_{j} against V_{j}. RESID contains OP*v_{j}. See STEP 3.
    call cgemv('N', n, j, -one, v, n, h(1, j), 1, one, resid, 1)

    if (j > 1) h(j, j - 1) = cmplx(betaj, 0.0)

    orth1 = .true.
    call ccopy(n, resid, 1, workd(ipj), 1)
70  continue
    ! Back from reverse communication if ORTH1 = .true. WORKD(IPJ:IPJ+N-1) := B*r_{j}.
    orth1 = .false.
    ! Compute the B-norm of r_{j}.
    rnorm = scnrm2(n, resid, 1)
    !-----------------------------------------------------------%
    ! STEP 5: Re-orthogonalization / Iterative refinement phase |
    ! Maximum NITER_ITREF tries.                                |
    !                                                           |
    !          s      = V_{j}^T * B * r_{j}                     |
    !          r_{j}  = r_{j} - V_{j}*s                         |
    !          alphaj = alphaj + s_{j}                          |
    !                                                           |
    ! The stopping criteria used for iterative refinement is    |
    ! discussed in Parlett's book SEP, page 107 and in Gragg &  |
    ! Reichel ACM TOMS paper; Algorithm 686, Dec. 1990.         |
    ! Determine if we need to correct the residual. The goal is |
    ! to enforce ||v(:,1:j)^T * r_{j}|| <= eps * || r_{j} ||  |
    ! The following test determines whether the sine of the     |
    ! angle between  OP*x and the computed residual is less     |
    ! than or equal to 0.717.                                   |
    !-----------------------------------------------------------%
    if (rnorm > 0.717 * wnorm) go to 100
    iter = 0
    !---------------------------------------------------%
    ! Enter the Iterative refinement phase. If further  |
    ! refinement is necessary, loop back here. The loop |
    ! variable is ITER. Perform a step of Classical     |
    ! Gram-Schmidt using all the Arnoldi vectors V_{j}  |
    !---------------------------------------------------%
80  continue
    ! Compute V_{j}^T * B * r_{j}.
    ! WORKD(IRJ:IRJ+J-1) = v(:,1:J)'*WORKD(IPJ:IPJ+N-1).
    call cgemv('C', n, j, one, v, n, workd(ipj), 1, zero, workd(irj), 1)
    ! Compute the correction to the residual:     |
    ! r_{j} = r_{j} - V_{j} * WORKD(IRJ:IRJ+J-1). |
    ! The correction to H is v(:,1:J)*H(1:J,1:J)  |
    ! + v(:,1:J)*WORKD(IRJ:IRJ+J-1)*e'_j.         |
    call cgemv('N', n, j, -one, v, n, workd(irj), 1, one, resid, 1)
    call caxpy(j, one, workd(irj), 1, h(1, j), 1)
    orth2 = .true.
    call ccopy (n, resid, 1, workd(ipj), 1)
90  continue
    ! Back from reverse communication if ORTH2 = .true.
    ! Compute the B-norm of the corrected residual r_{j}.
    rnorm1 = scnrm2(n, resid, 1)
    ! Determine if we need to perform another step of re-orthogonalization.
    if (rnorm1 > 0.717 * rnorm) then
        !---------------------------------------%
        ! No need for further refinement.       |
        ! The cosine of the angle between the   |
        ! corrected residual vector and the old |
        ! residual vector is greater than 0.717 |
        ! In other words the corrected residual |
        ! and the old residual vector share an  |
        ! angle of less than arcCOS(0.717)      |
        !---------------------------------------%
        rnorm = rnorm1
    else
        ! Another step of iterative refinement step is required. NITREF is used by stat.h
        rnorm = rnorm1
        iter = iter + 1
        if (iter <= 1) go to 80
        ! Otherwise RESID is numerically in the span of V
        do jj = 1, n
            resid(jj) = 0.0
        end do
        rnorm = 0.0
    end if
    ! Branch here directly if iterative refinement
    ! wasn't necessary or after at most NITER_REF
    ! steps of iterative refinement.
100 continue
    rstart = .false.
    orth2 = .false.
    ! STEP 6: Update  j = j+1;  Continue
    j = j + 1
    if (j > np) then
        ido = 99
        do i = max(1, k), np - 1
            ! Check for splitting and deflation.
            ! Use a standard test as in the QR algorithm
            ! REFERENCE: LAPACK subroutine zlahqr
            tst1 = slapy2(real(h(i, i)), aimag(h(i, i))) + slapy2(real(h(i + 1, i + 1)), aimag(h(i + 1, i + 1)))
            if(tst1==real(zero)) tst1 = clanhs('1', np, h, ldh, workd(n + 1))
            if(slapy2(real(h(i + 1, i)), aimag(h(i + 1, i))) <= max(ulp * tst1, smlnum)) h(i + 1, i) = 0.0
        end do
        return
    end if
    ! Loop back to extend the factorization by another step.
    go to 1000
end
